Lattice operators from discrete hydrodynamics 
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We present a general scheme to derive lattice differential operators from the discrete velocities 
and associated Maxwell-Boltzmann distributions used in lattice hydrodynamics. Such discretizations 
offer built-in isotropy and recursive techniques to increase the convergence order. This provides a 
simple and elegant procedure to derive isotropic and accurate discretizations of differential operators, 
which are expected to apply across a broad range of problems in computational physics. 



The development of efficient and accurate lattice ver- 
sions of differential operators is a central theme of compu- 
tational physics. Indeed, the numerical solution of any 
classical or quantum field theory requires the develop- 
ment of discrete differential operators in order to solve 
the partial differential equations associated with the con- 
tinuum theory. Isotropy of the differential operators in 
continuum space is often lost in developing the corre- 
sponding discrete operators. This inability to perform 
discrete operations satisfying the inherent isotropy of 
continuum space may reflect severely on numerical sim- 
ulations of physical problems. It is therefore desirable to 
develop discrete operators which retain as many symme- 
tries as possible of their continuum counterparts. Here, 
we address this issue with specific regard to isotropy, and 
show that progress in lattice hydrodynamic simulations, 
naturally provides a strategy to develop such operators 
in full generality, beyond the original realm of hydrody- 
namics. 

Finite difference schemes remain one of the most popu- 
lar method of discretising differential operators. Though 
the accuracy of the scheme can be improved by increas- 
ing the stencil size, discrete operations are generally re- 
strained to the principal directions of the lattice (coor- 
dinate directions on a rectangular grid), often neglecting 
the grid points along other directions. For example, in or- 
der to calculate the curl or Laplacian of a field, very often, 
only information available at principal directions is used. 
This leads to a loss of information which deteriorates the 
accuracy of the discrete operation, isotropy in the first 
place. Use of larger stencils with next-nearest neighbors 
may offer significant improvements without degrading ef- 
ficiency to any significant extent. This issue has been 
addressed before in the literature. For instance, mimetic 
discretizations have been developed to recover the prop- 
erties of underlying continuum theory [1]. The isotropy 
of Laplacian operators has been the object of several pre- 
vious studies [2] and a specific illustration of the general 
method proposed was presented in [3]. In this Letter, we 
show that the same procedure can be applied in full gen- 
erality to a broad class of differential operators, such as 
gradient, divergence and curl, which play a central role 
across virtually all of computational physics. 



Discrete operators from lattice kinetic theory: Consider 
a unit cell of dimensions 2Ax x 2Ax x 2Ax as shown in 
Fig. (1), generating a standard uniform grid in Cartesian 
coordinates. The center of the cube is the point of in- 
terest and will be denoted as r. Neighboring points on 
the grid vary in distance and can be classified as near- 
est neighbors - NN, next nearest neighbors - NNN, and 
next-next-nearest neighbors - NNNN, as highlighted in 
Fig. (1). Let the vectors pointing to each of these points 
be denoted by c^, where i represents any point on the 
grid. We include Co which is a zero vector at r. Let 
us also define weights, Wi, associated with these different 
points on the lattice. To be more precise, within the con- 
text of lattice kinetic theory, represent a set of discrete 
speeds which move the information from the center point 
r to the z-th neighbor r + in a time-step At, according 
to the light-cone rule, = c^Ai. In the following, we 
shall take At = 1, so that the discrete displacements 
can be identified with the discrete speeds Cj. 

We begin with a DnQm lattice hydrodynamic model in 
n dimensions with m discrete velocities. It is well known 
from the lattice hydrodynamics literature that, in order 
to preserve isotropy of discrete space up to fourth order 
in the lattice tensors, it is necessary to have [4] 



J2 Wi = 1 

i 



^WiCi^ a Ci^c in Ci t x = T 2 A {i ' 



(1) 
(2) 
(3) 



where Greek indices label Cartesian directions and 

^al-yX = fW^A + 5 a \S~fP + S a7 5/3\. We cllOOSC Ax = 1 

such that T = 1/3, a lattice-dependent constant. It is 
well known that the lattice formulation of kinetic theory 
provides a computationally efficient method to solve con- 
servation equations and is well established for the hydro- 
dynamics [5] . A discrete form of the Maxwell-Boltzmann 
velocity distribution [6] is used in this formulation, which 
preserves the isotropy of space to fourth order. Higher 
order stencils may be used to obtain 6th order accuracy 
[7]- 

We now introduce the method of generating various 
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FIG. 1. Points on a cubic cell of dimensions 2Aa; x 2Ax x 2Ax. 
is the point of interest and will be represented by r. Here NN 
represent the nearest neighbors, NNN the next nearest neigh- 
bors, and NNNN the next next nearest neighbors, marked by 
o, •, x and ■ respectively. 



discrete operators which preserves isotropy upto fourth 
order. Isotropy at higher orders needs larger stencils. 
Now consider the general operator (s = V) 



■/s - J 



(4) 



that acts on a field ip(r) discretely defined on a lattice as 
shown in Fig. (1). As a result, 

fs Mr)] = i E me^W = w ^ r + c ^ ^ 



We illustrate below how / s can act as a generating func- 
tion to construct several discrete operators. 

Consider the following two transformations = 



2 (f a - fo) ip(r) and Q(tp) = %ip(r) or equivalcntly 



N 



(6) 



1=1 

N 



G($) = if, E WiCii/j(r + a) 



where N is the total number of neighboring points con- 
sidered. Taylor expanding these expressions and using 
the symmetries of Eqs. (1 - 3), we obtain 

J"(V>) = ( 1 + J v2 ) V ^ and = ^1 + ^V 2 ^) Vtp. 

These expressions may be solved for V 2 </> and V</> re- 
spectively by inverting the linear operators. We retain 
only the leading order terms and show how these expres- 
sions may be used to obtain isotropic discrete operators 
and how a recursive technique can be employed to obtain 
higher order accurate discrete operators. 



Considering only the leading order terms of the in- 
verted linear operators , we find that 



V^V = - 

v T 



N 



J~]wii/)(T + Cj) -ip(r) 



N 



0(V 4 



V ^ = f £ mc t ifj(r + a) + 0(V 3 ). 



(8) 



(9) 



These expressions preserve isotropy at the leading or- 
der. The former expression has been reported earlier in 
[3]. The latter is also suggested and used in various con- 
texts [8]. Here we explain the underlying connection to 
lattice hydrodynamics and present the corresponding ex- 
pressions for the divergence and the curl. Hence, for a 
vector field, <j>(r), 



V-0=i£ W4 c 4 -</)(r + c J ) + O(V 3 ) (10) 

i=l 
1 N 

V A = - 2J WiCi A </>(r + Ct) + 0(V 3 ). (11) 

i=l 

These expressions are remarkable because they display 
isotropy up to leading order error, with an error coeffi- 
cient of order 0(T). Any lattice with suitable weights 
which satisfies the conditions in Eq. (1-3) will provide 
an expression for the discrete operators and will ensure 
isotropy. For example the weights used in the lattice 
hydrodynamics may be used. These weights are lattice 
analogues of the Maxwell-Boltzmann equilibrium in con- 
tinuum velocity space. 

A recursive algorithm can be developed now to remark- 
ably improve the accuracy of these operators. Consider- 
ing only the next order terms of the inverted linear op- 
erators , we have 



(7) In the explicit form, these expressions are 



(12) 
(13) 
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T 
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y Wjip(r + Cj) - *0(r) 

i=l 
N 

^^Wr + c,))-^)) 



(14) 



1 N 

= -^] WiCi4'(r + Cj) 



i=l 
N 



J2^Wr + Cl ))~g^(r)) 



i=l 



(15) 



By employing this approximation for the next order term, 
we are essentially following a predictor-corrector method, 
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N, (for 2D) 


D2Q9 


D3Q15 


D3Q19 


D3Q27 





1 (1) 


4/9 


2/9 


1/3 


8/27 


NN 


6 (4) 


1/9 


1/9 


1/18 


2/27 


NNN 


12 (4) 


1/36 





1/36 


1/54 


NNNN 


8 (0) 





1/72 





1/216 



TABLE I. Popular models and the associated weight factors 
for various DnQm lattice hydrodynamics models. Values of 
N for the two-dimensional D2Q9 model are given in brackets. 



however with a much simpler implementation. The trans- 
formation Q(ij>) can be easily manipulated to obtain the 
divergence and curl operators. For a vector function </>(r) 
defined on a grid, the generating function can be manip- 
ulated as V(<j>) = 4- ■ 0( r ) and C (0) = A <M r ) to 
as as 

obtain 



ds 



V(<j>) = ^J2w t c l -cf>(r + c l ) (16) 

i=i 

1 N 

C(<t>) = 7f;Ys WiCi A0(r + Ci). (17) 



which provide the expressions for divergence and curl: 



v • = dm - 

VA0 = C(V)-|^(C(V)). 



(18) 
(19) 



By employing the recursive technique we obtain a higher 
order accurate method but the leading order error, which 
is 0(V 5 ) in case of Eq. (15), (18) and (19), is not 
isotropic. As mentioned earlier, it is not possible to get 
the isotropic error beyond fourth order with this lattice. 
However, it is possible to go to larger stencils and ob- 
tain isotropy at the desired level in the same formulation. 
Double-differential operators, such as V(V-), V • (VA), 
V A (VA) and similar, may be derived likewise from the 
generating function using d 2 //ds 2 . 

In short, for any transformation K,(4>) which at leading 
order provides a differential operation on <f>, say if [</>], it 
is possible to write 



K[ct>)^K{ct>)-XF[lC{<j>)] 



(20) 



providing a method to perform the same differential op- 
eration. In the above expression, A = T/4 for the Lapla- 
cian, and A = T/2 for gradient, divergence and curl op- 
erators. Discrete operators can now be derived using any 
standard DnQm models. Popular models and the as- 
sociated weights in lattice hydrodynamics literature are 
listed in Table (I). 

Results and Discussion: We next compare the ac- 
curacy and isotropy properties of the discrete op- 
erators, Eqns. (14 - 15, 18 - 19) that were de- 
rived in the last section. The discrete Fourier trans- 
form of the gradient operator is given by D(k) = 



J2 r exp(— ik.r)Q (ip) / ^ r exp(— ik.r)^(r). In the small 
wavelength limit, the corresponding expressions for dif- 
ferent models and the standard second order central dif- 
ference (CD) scheme may be written as follows. For clar- 
ity, only one component is shown below. 

1 _______ + (k 6 ) 

36 180 1 ' 




where k = |k| and (a ^ (3 =/= 7) are the cartesian in- 
dices. For curl and divergence, the form of these oper- 
ators remains the same. Note that repeated indices are 
not summed upon. 

Contour plots of these operators in Fourier space are 
plotted in Fig. (2). The error involved in these calcula- 
tions may be estimated by comparing with the analytical 
value |k|. Our schemes always provide better isotropy 
compared to the standard finite difference method at 
small wave numbers. Such a comparison for Laplacian 
operators can be found in [3]. We now apply these de- 
rived discrete differential operators to various test func- 
tions and compare them with standard difference meth- 
ods in the literature. 

We consider a two dimensional Gaussian as a test func- 
tion, 



ip(x, y) = exp 



-(x 2 +y 2 ) 



2a 2 



(21) 



where a is the variance. In a square domain of length 
unity spanned by N g grid points in each direction, 
we compute the gradient using Eq. (9) for a D2Q9 
model. The gradient is also calculated using the mod- 
ified scheme, Eq. (15), and the standard second order 
central difference formula. Comparing with the analyti- 
cal expression for the gradient, V^>(x, y) anal y t ) L 2 — norm 
of the error is defined as 



norm 



v Ef [v W D2 " - vw 



analytl 



(22) 



The results are illustrated in Fig. 3(a), where the error is 
plotted as a function of grid spacing, suitably nondimen- 
sionalised. The increase in accuracy with the increase 
in number of grid points is apparent. The second order 
convergence of Eq. (9) and the fourth order convergence 
of Eq. (15) may also be seen. While our lower order 
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(d) D3Q19 (e) D3Q27 

FIG. 2. Contour plots in Fourier space of the discrete gradient operator D(k), are shown for the planes k z = and k z = 1 for 
(a) continuous (b) central difference, (c) D3Q15, (d) D3Q19 and (e) D3Q27 models. The present schemes invariably exhibit 
better isotropy. For clarity, the y-axis is omitted in k z = 1 plots. 
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(a) Gradient (b) Curl (c) Divergence 



FIG. 3. L 2 -norms of the error for (a) discrete gradient operator applied on a Gaussian function, Eq. (21), (b) curl operator on 
Taylor-Green velocity field, Eq. (23), and (c) divergence operator on a Gaussian vector field, in two dimensions are plotted as 
a function of grid spacing in log-log scale. Since the domain is of fixed size of unity for (a) and (c), abscissa is suitably scaled 
as a/ Ax = aN g . We have used a = 0.05. In case of (b), k x = 3, k y = 4 and domain length = 2n are chosen, hence, abscissa is 
nondimensionalised as (2n/k)/Ax. For clarity, the ordinate is labelled only in the left most panel. Dashed lines of slope 2 and 
4 are plotted for eye-guiding purposes. The significant improvement of the recursive scheme may be noted in all cases. 



scheme maintains the isotropy property compared to the 
standard central difference method, the recursive scheme 
improves the accuracy considerably, i.e to fourth order. 

Next we address the curl operator, whose discretiza- 
tion is central not only to hydrodynamics but also to 
computational electromagnetics [9]. As a test function, 
u(x,y), we choose, 



u x = -k y cos(k x x) sm(k y y) 
u y = k x cos(k y y) sin(k x x) 




which represents the velocity field in a Taylor-Green flow. 
The vorticity, which is the curl of velocity field, is given 
by fc 2 cos(k x x) cos(k y y)z where fc 2 = fc 2 + ky. The curl 
of this field is also obtained for a D2Q9 lattice model, 



using both lower order scheme and the improved scheme, 
Eqs. (11) and (19) and compared with the analytical ex- 
pression. The error, defined using L 2 — norm as earlier, 
is plotted in Fig. 3(b). The curl, as calculated using the 
standard central difference scheme, is also plotted in the 
same figure. It may be noted that, like for the gradient, 
the discrete expressions derived here provide significantly 
higher accuracy. Again, it may be noted that isotropy is 
in-built in the lower order operator. 

We now compute the divergence of a Gaussian vector 
field, Vf/;(a;,?/) analyt using the D2Q9 model. Equations 
(10) and (18) were used to compute the divergence of 
this function and compared with the analytical expres- 
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sion exp [— (x + y )/2er ] j ^ . The corre- 
sponding error, as L 2 norm, is plotted in Fig. 3(c). The 
improvement in accuracy, as compared to a standard cen- 
tral difference operator, is also illustrated. 

It may be mentioned that the Laplacian of a Gaus- 
sian function is calculated and the error, as Li norm, for 
the discrete operators, Eqns. (8) and (14), behaved in a 
similar fashion, witnessing the generality of the method. 

Summarizing, we have shown that the stencils as- 
sociated with the discrete-velocity schemes developed 
in lattice hydrodynamics, naturally provide an elegant, 
efficient and accurate procedure to formulate discrete 
isotropic versions of the most fundamental differential 
operators, such as gradient, curl, divergence and Lapla- 
cian. Furthermore, we have also shown that the accu- 
racy of these operators can be systematically improved 
by means of a recursive iteration procedure. Applica- 
tion of these discrete operators to various smooth test 
functions, was shown to result in significantly improved 
accuracy, as compared to standard finite-difference op- 
erators. Finally, we wish to emphasize that, owing to 
its generality, the present method is expected to apply 
to further classes of differential operators, such as the 
Dirac propagator and Wilson plaqucttcs in lattice gauge 
theories [10]. 
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